import numpy as np
from scipy.spatial.distance import pdist, squareform
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs

# 生成测试数据
def generate_test_data():
    # 生成三个聚类
    centers = [[0, 0], [3, 3], [6, 0]]
    X, y = make_blobs(n_samples=1300, centers=centers, cluster_std=0.5, random_state=42)
    
    # 生成噪声点（30个）
    x_min, x_max = X[:, 0].min()-1, X[:, 0].max()+1
    y_min, y_max = X[:, 1].min()-1, X[:, 1].max()+1
    
    noise_x = np.random.uniform(x_min, x_max, size=30)
    noise_y = np.random.uniform(y_min, y_max, size=30)
    noise = np.column_stack((noise_x, noise_y))
    
    # 合并数据
    data = np.vstack([X, noise])
    true_labels = np.hstack([y, np.full(30, -1)])  # 噪声标签为-1
    
    return data, true_labels

# 密度峰值聚类实现
def density_peak_clustering(data, n_clusters=3, dc_percent=2.0, noise_percent=5):
    # 计算距离矩阵
    distances = squareform(pdist(data))
    n = data.shape[0]
    
    # 自动确定截断距离dc
    num_pairs = n * (n - 1) // 2
    dist_list = np.sort(distances[np.triu_indices(n, k=1)])
    dc = dist_list[int(dc_percent/100 * num_pairs)]
    
    # 计算局部密度rho
    rho = np.sum(np.exp(-(distances / dc)**2), axis=1) - 1
    
    # 计算距离delta和最近更高密度点
    order = np.argsort(-rho)
    delta = np.zeros(n)
    nearest_higher = np.zeros(n, dtype=int)
    
    # 处理最高密度点
    delta[order[0]] = np.max(distances[order[0]])
    nearest_higher[order[0]] = -1
    
    for i in range(1, n):
        current_idx = order[i]
        higher_rho = order[:i]
        dists = distances[current_idx, higher_rho]
        if len(dists) > 0:
            min_idx = np.argmin(dists)
            delta[current_idx] = dists[min_idx]
            nearest_higher[current_idx] = higher_rho[min_idx]
        else:
            delta[current_idx] = 0.0
            nearest_higher[current_idx] = -1
    
    # 选择聚类中心
    product = rho * delta
    centers = np.argsort(-product)[:n_clusters]
    
    # 分配类别
    labels = np.full(n, -1, dtype=int)
    for cluster_id, center in enumerate(centers):
        labels[center] = cluster_id
    
    for idx in order:
        if labels[idx] == -1:
            labels[idx] = labels[nearest_higher[idx]]  # 补全赋值语句
    
    # 标记低密度点为噪声
    if noise_percent > 0:
        rho_threshold = np.percentile(rho, noise_percent)
        labels[rho < rho_threshold] = -1
    
    return labels, centers, rho, delta

# 生成测试数据
data, true_labels = generate_test_data()

# 执行聚类
labels, centers, rho, delta = density_peak_clustering(data, noise_percent=5)

# 可视化结果
plt.figure(figsize=(15, 6))

# 真实分布
plt.subplot(1, 3, 1)
plt.scatter(data[:, 0], data[:, 1], c=true_labels, cmap='viridis', s=10)
plt.title("True Distribution with Noise")

# 决策图
plt.subplot(1, 3, 2)
plt.scatter(rho, delta, s=10)
plt.scatter(rho[centers], delta[centers], c='red', s=50, marker='X')
plt.xlabel('Density (rho)')
plt.ylabel('Distance (delta)')
plt.title("Decision Graph")

# 聚类结果
plt.subplot(1, 3, 3)
plt.scatter(data[:, 0], data[:, 1], c=labels, cmap='viridis', s=10)
plt.scatter(data[centers, 0], data[centers, 1], c='red', marker='X', s=100)
plt.title("Clustering Result")

plt.tight_layout()
plt.show()
